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Abstract 

The generating functional method is employed to investigate the synchronous dynamics of 
Boolean networks, providing an exact result for the system dynamics via a set of macroscopic 
order parameters. The topology of the networks studied and its constituent Boolean functions rep- 
resent the system's quenched disorder and are sampled from a given distribution. The framework 
accommodates a variety of topologies and Boolean function distributions and can be used to study 
both the noisy and noiseless regimes; it enables one to calculate correlation functions at different 
times that are inaccessible via commonly used approximations. It is also used to determine condi- 
tions for the annealed approximation to be valid, explore phases of the system under different levels 
of noise and obtain results for models with strong memory effects, where existing approximations 
break down. Links between BN and general Boolean formulas are identified and common results 
to both system types are highlighted. 
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I. INTRODUCTION 



Boolean networks (BN) have been suggested as simplified models of various biologieal 

systems, in particular for modeling gene- regulatory networks Simplifying the state of 
genes by adopting a two state variable representation, ON/OFF, enables one to model the 
complex interactions between genes as arbitrary Boolean functions of randomly selected 
variables. This model gives rise to a rich and complex behavior that has been successfully 
employed to gain insight into the dynamics and steady states of various gene-regulatory 
systems 2I. 

BN models comprise sites (genes), each of which is represented by a binary variable. 
The state of each variable is determined by the states of k randomly selected sites via a 
fc-input Boolean function. The specific k input variables selected for each site constitutes 
the network topology and the selected Boolean functions determine the corresponding in- 
teraction leading to the variables' state. In the original formulation |]]], both topology and 
Boolean functions were selected uniformly from the ensemble of networks with in-degree 
k per variable and 2^*" possible Boolean functions, respectively. Both are considered fixed 
{quenched variables). Moreover, the original model was deterministic and the analysis there- 
fore focused on its periodic-orbit attractors, steady-state and their basins of attraction. 

While this family of networks was originally introduced to model the gene-regulatory 
network and is commonly known as Random Boolean Networks (RBN) or Kauffman 
nets, similar topologies have been employed to study network properties in other application 



domains, ranging from social [3] to genetic and neural [5| networks. Although the topology 
used is common to all these models, based on a discrete state-space and random fc-variable 
Boolean interactions, the nature of the interaction may be different for each of the models. 
We will refer to the general class of A^-variable binary system with connectivity degree k as 
the N-k model [g,^. 

This abstraction of complex gene-regulatory system lends itself to analysis in terms of 
both their dynamics and equilibrium properties ja, S]. Equilibrium analysis relies mainly 
on the cavity method, while the dynamics has been mostly investigated using the annealed 
_„«„„Qduetoitss,™pl,eitvaud success in providing accurate results in many of 
the models studied, especially for very large systems. The underlying approximation in this 
approach is that both thermal and quenched variables (primarily the network topology) are 
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considered to be on equal footing and are sampled at each time step. This helps suppress 
emerging correlations of specific sites at different times, simplifies the analysis and gives rise 
to an effective methodology which works in most cases. The annealed approximation was 
particularly successful in large-scale systems (A^ — )■ oo); it allows one to predict the evolu- 
tion of network activity and Hamming distance order parameters. The former refers to the 
magnetization or the proportion of ON/OFF states and the latter to the difference between 



the states of the network starting from different initial conditions. It was shown [11|, [12 1 
that the annealed approximation provides accurate magnetization and Hamming distance 
order parameter predictions for RBN (i.e., with uniformly sampled Boolean functions); how- 
ever, the conditions for its applicability and validity for general N-k systems has remained 
unclear ^|. Moreover, in some cases, especially in systems with memory, discrepancies 
have been found between results obtained via the annealed approximation and simulation 
results [l^, casting doubts on the validity of the approach for such models. 

The aim of the current paper is to develop further a framework for exact analysis of 



N-k models based on the generating functional analysis (GFA) framework 15|, |l6|, which 



has been employed successfully in the study of various Ising spin models 17|. The newly 
developed framework is then employed to determine the conditions under which the annealed 
approximation is valid, to investigate the possible phases of BN depending on the noise level 
and to demonstrate its efficacy for analyzing systems with memory, where the annealed 
approximation is known to break down [14]. We note that an alternative to the GFA 
method (called the dynamical cavity method) was recently introduced 18| , which we believe 
can be also used for the range of models studied here. 

Section Hi] introduces the BN model while Section HH] describes the methodology used 
and its application to the current model. In section HVl we employ the dynamical equations 
obtained for the order parameters to investigate the conditions under which the annealed 
approximation provides exact results; a similar set of equations is then used to identify 
possible phases of the system and to analyze the dynamics of system with memory. Finally, 
in section |Vl we summarize the results obtained and point to future research directions. 
Some of the detailed derivations appear in dedicated appendices. 



3 



II. MODEL 



The model we consider here is a recurrent Boolean network which consists of binary 
variables Si(t) G {—1, 1} interacting via Boolean functions : {—1, 1}^ — )• {—1, 1} of exactly 
k inputs. Because of thermal noise, which can flip the output of a function with probability 
p (l9|, a site Si(t) in the network is operating according to the stochastic rule 

S,{t + l) = r],{t)a,{Si,{t),...,S,,{t)), (1) 

where r]i{t) is an independent random variable from the distribution P{r]) = + — 

p)S>ri-i. The function-output ^^(t + l) is completely random when p = 1/2 and completely 
deterministic when p = 0. Averaging out the thermal noise {f]i{t)} in the system governed 
by ([1]) gives rise to the microscopic law 

l3Si{t+l)ai{Si^{t),..,Si^{t)) 

where the inverse temperature /3 = 1/T relates to the noise parameter p via tanh/3 = l — 2p. 

All sites in the network are updated in parallel and given the state of the network S(t) G 
{ — 1,1}^ at time t, the function-outputs ^(t + l) at time t + 1 for the different sites are 
independent of each other. This Markovian property allows us to write the probability of 
the microscopic path S{0) —)■■■■—)■ S{tmax) as a product of over all sites and time 
steps. Furthermore, we consider two copies of the same topology but with different initial 
conditions, shown in Figure [H comparing the two will enable us to study the effects of 
initial-state perturbations. Following similar arguments to those of the single network case, 
the joint probability of microscopic states in the two systems are given by 

P[{5(t)};{5(t)}] = P(5(0),5(0)) J] P(5(t+l)|5(t))P(5(t+l)|5(t)), (3) 

t=o 

where P{S{t+l)\S{t)) = ]\l,P^XS^it + l)\S,At). --.S^M)- 

The sources of quenched disorder in our model are random Boolean functions and random 
connections. Boolean functions {at} are sampled randomly and independently from the 
distribution 

P{(y) = ^Pi K,o^, (4) 
7eB 
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FIG. 1: The model of two systems with identical topology and different initial conditions. Functions 
are indicated by squares and input nodes by circles. White (striped) indicates flipped inputs. 

where Yl-yesPi = 1; P7 ^ and B is the set of all fc-ary Boolean functions. The connectivity 
disorder arises from the random sampling of connections generated by selecting the i-th 
function and sampling exactly k indices, I = {ii, .., ik}, uniformly from the set of all possible 
indices [A^] = {!,..., N}. This gives rise to the probabilities 
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where Za is a normalization constant. The connectivity tensors {A-^ j^} define the random 
topology via entering into the definition of probability with ai{Si-^(t), .., Si^^it)) being 
replaced by i^. ••, Si^(t)). Other connectivity and function profiles can 

be easily accommodated within our framework by incorporating additional constraints into 
the definitions (jl]) and ([5]) via the appropriate delta functions. 

III. GENERATING FUNCTIONAL ANALYSIS 



To analyze the typical properties o f t 
functional method of De Dominicis 
the generating function 



le system governed by ([T]) we will use the generating 



15j. Following the prescription of 15| we first define 



(6) 



where (...) denotes the average over all paths occurring in two systems governed by the 
joint probability (E]). The generating function ([2D allows us to compute moments of 
([3]) by taking partial derivatives with respect to the generating fields {'ijji{t),ipj{s)}, e.g. 



{Si{t)Sj{s)) = — lim^^ (^^d- , — -'^f''/'' V']- Secondly, we assume that the system be- 
comes self-averaging, i.e. all thermodynamic macroscopic properties are self-averaging, for 
— 7- oo [15] and compute r['0; -0], where is the disorder average; this gives rise to the 
macroscopic observables 



w4i:Wi)> = n.n (7) 



where m{t) is the network activity (or magnetization (20I), C(t,s) is the correlation between 
two states of the same network and C^it) is the overlap between two copies of the same 
network which is related to the Hamming distance d{t) via d(t) = — Ci2(t)). 

Averaging the generating function (|6]) over the disorder |35l] (see Appendix |A] for details) 
leads us to the saddle-point integral problem 

T = j {dPdPd(]dn}e^*[^^'^'^'^>l (8) 

where \1/ is the macroscopic saddle-point surface 

^ = ij2 pi{sm)p{{s{t)}) (9) 



+i J {dh{t)} duj n{{h{t)},uj) n{{h{t)},uj) 
+ [{dh{t)}doon{{h{t)},uj) J2 IflPiiSAt)}) 

^ {S,{t)} li=l 

+ log J2 /{dh(t)dh(t)}^'^^M[{5(t)},{h(t)}|{h(t)},a;,{0}] 

using the notation S{t) = {S^{t), S^{t)), h{t) = {h\t) , h"^ (t)) , h{t) = {h\t),h'^{t)); M is an 
effective single-site measure 



M[{5(t)}, {h(t)}|{h(t)}, CO, {0}] = P{S\0), S'm e 



-iQ{{h{t)},Lu)+iiu-iP{{S (t)}) 



tmax-1 2 ^ „/3ST(i+l)/l7(i) 



with P{S\0), ^2(0)) = 1(1 + S\0)m{0) + S^{0)m(0) + ^^(0)^2^0)^12(0)). The generating 
fields 1/), have been removed from the above as they are not needed in the remainder of this 
calculation. In the limit of — t- 00 the integral ([8]) is dominated by the extremum points 
of the functional \1'. The functional variation of \E' with respect to the order parameters 
{P, P, Q, Q} leads us to the saddle-point equations 

/ tmax 1 \ 

P{{Sit)}) = i n (11) 

k ( k \ 

p{{s{t)}) = iY, E ^[{sm{s^m]\\[p{{s,m\ (12) 

*=l{S,(t)} Ij^i J 



X 



j {dh(t)} da; n{{h{t)},u) e-'^'-o-''J:Uf^'(i) -is'lit),-,s2(t))-i. 



nam}, CO) 



t=0 
k 



6{uj - uj') ) (13) 

M 



Q{{m},uj) = iY, J]P({5,(t)})e-'»^^^-'^^W"(^iW'-'^^'«)-^- , (14) 
{s,(t)}i=i 

where (. . is average generated from the single-site measure flTUl) . In Appendix [B] we 
show that the conjugate order parameter P is a constant. Using this result the saddle-point 
equation ([Hj) in the single-site measure (fTO|) leads us to the main equation of this paper 

k 

PiS,S) = P{S{0),SiO))J2 l[[P{s„s,) 



i-max 1 

X 

i=0 



l[P^iSit+l)\Siit), .., Skit)) P^iSit+l)\Siit), .., Skit)) . (15) 



The physical meaning of ( ITS!) is revealed by PiS, S) = limjv->oo jf Xlili i^'^i •^i] 

the disorder- averaged joint probability of single-spin trajectories S and S in the two systems. 

Equation (HM can be used to compute the macroscopic observables ([7]). To demonstrate how 
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this can be done we derive explicitly the expression for the two time correlation C{t, s) 



5^P(5,5) S{t'+1) S{t"+1) 



(16) 



s.s 



J]p(5(o),^(o)) llns,,s,)\s{t'+i)s{t"+i) 

{s„s^}i=i 



s,s 



tmax 1 



X II P„(5(t+l)|5i(t),..,S,(t))P„(5(t+l)|Si(t),..,^fe(t)) 
i=0 

k 

{Sj{t'), Sj{t"}j=l 

X PaiSit'+l)\S,{t'), .., Skit')) P^{S{t"+l)\S,{t"), .., " 

S{t'+l),S{t"+l) 

X 5(t'+l) 5(t"+l) = ; 

note that many of the variables in the summation over P{S, S) are redundant, they have 
been introduced for methodological reasons but vanish during the derivation. Carrying out 
a similar derivation for the other order parameters one obtains a closed set of iterative 
equations: 

k 



m 



(t+l) = /„(m(t)) = tanh(/3) 

s j=l 

C{t + 1, s + = P„(m(t), m{s), C{t, s)) 



1 + Sjm{t) 



a{S) 



= tanh2(/3)^n 

Ci2(t + 1) = P„(m(t),m(t),Ci2(t)), 



l + Sjm(t) + Sjm{s)+SjSjC{t, s) 



a{S) a{S) 



(17) 

(18) 
(19) 



where S = {Si, . . . , 5*^) and the equation for rh(t) is the same as (fTTI) . 



IV. RESULTS 

In this section, we first apply the equations f|T7|) - f|T9|) to the recurrent Boolean networks 
with thermal noise. We recover results of the annealed approximation for the order param- 
eters m and Ci2- However, the two-time correlation function C{t, s), computed here for the 
first time, allows us to study properties of the stationary states. Furthermore, the exactness 
of our method allows us to derive a rigorous upper bound on the noise level above which 
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FIG. 2: Topology of the basin of attraction with its stationary states. 

the system is always ergodic. In addition, we use the equation (1T5|) to study models with 
strong memory effects where the annealed approximation method is no longer valid. 

A. Boolean networks 

1. Stationary states 

It is clear from the results ffT7|l - ffT^ and f[T^ that the evolution of all many-time single-site 
correlation functions is driven by the magnetization m{t). A similar scenario was observed in 
recurrent asymmetric neural networks [17] which have a similar topology but uses different 
update functions than model ([T]) . The asymmetric neural network model can be seen as 
a special case of the N-k model when only linear threshold Boolean functions are used and 
the thermal noise enters into the system via randomness in the thresholds. Furthermore, 
for the stationary solution m = (m = limj__j.oo '^(^)) the solution of g = Fa{m,m,q) 

[here q = \imt^oo^i^,t„^-^ooC{t + tuj,tuj) is the Edwards- Anderson order parameter, used in 
disordered systems [21] to detect the spin glass (SG) phase where m = and q 0] and 
Ci2 = Fa{m, m, C12) are identical. This was also observed in asymmetric neural networks 



22 



and because of the equality q = Cu there is only one average distance ^(1 — g) on the 



attractor 



22| and all points in the basin of attraction uniformly cover the stationary states 



(see Figure [2]). 
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2. Annealed model 



The annealed approximation method, where connectivities and Boolean functions change 



at each time step of the process ([T]), provides identical results for m and Cu to those of ([Tj 
and f|T9l) [3]. However, within annealed approximation the two-time correlations take the 
form C{t, s) =m{t)m{s), where t>s, which is the solution of ( !T8|) only when networks are 
constructed from a single Boolean function. This result follows from the equality C{t, 0) = 
m{t)m{0), which is clear from the equations f llSp and fll6p . and the fact that for a single 
function, in the absence of an average over a, the joint probability of two spins in the 
equation f lTSj) factorizes when C{t, s) = m(t)m(s). 

The classical annealed approximation result 10 1 for RBN, which is exact in this case [ll[ 



12| . can be easily recovered from the equations f lT7|) - f|T9|) using the property a{S) =0 for 

~ — 

all S G { — 1, 1}''' and a{S)a{S) = 0, for all S S where the a average is taken over all 
Boolean functions with equal weight. In the noisy case (/3<oo), the magnetization variable 
m{t) = for all t > and q = tanh^(/3)(^)*'', corresponding to the stationary solution of 
( 1T8|) . has one stable solution q^O for all finite f3 > and k. For /3— j-oo (no noise), there is 
a transition from one stable solution q = 1 for k <2 to two solutions q = l (unstable) and 
q^O (stable) for fc>2 0. 

3. Noise upper bound 

An interesting question related to the ergodicity and phase transitions is whether the 
system ([T]) can retain information about its initial state in the presence of noise. This 



question has received a considerable attention in the works on cellular automata 
a closely related field of fault-tolerant computation 24 



23| and 



m 



The unordered paramagnetic (PM) phase m = 0, where no information is retained, is a 



fixed point of f lT7|) only when ^5.a;(S') =0. 

Proposition IV. 1 The point m = is a stable and unique solution of [Tl\ ) when tanh/3< 
2''"V^((fe-r}/2)!2''"V(^-l)((fe!2V2)}=K^) for k odd and even respectively. 

Proof To prove this we first find a Boolean function x such that fxim) > faim) when 
m G [0,1) and /^(m) > f^im) when m G (—1,0]. It turns out that any function from 
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the set x{S) = sgn[y^^ , Sj ] + S\0; T^^-i 5'j]7('5'), where 7(5') G {—1,1} and such that 
Ss'^iO; Sf=i =0 36| satisfies these properties. To show this we define the average 



\' ■ ■ )s\m ~ '^S T\.j=l 



1+Sim 



and use the shorthand notations {l+fS], l-i^], lo['S']} for 



the indicator functions Sj > 0], lEj=i Sj < 0], l[Xlj=i Sj = 0]}. Then we compute 

the difference A(m) = (/^(m) — /Q(m))/4tanh(/3) as foUows 



.j=i 



-a{S) 



(20) 



Sim 



4(1+ 



:i+ [5] + l_[5] + lo [S])a{S) 



S\m 



-^1+[S]-(1-«(S) )-l_[S]-{l+a{S) )--l,[S]a{S) 

■a 1 



-(U[S]l[aiS) = -l] -1_[5]1[«(5) = +1] --lo[S]aiS) 



S\m 



S\m 



-E 



1 + m 



l — m 
1+m 



l — m 



1+m 
l—m 



U [S]l[a{S) = -l] 



l_[5]lK5) = +l]"-i5^1o[5]a(5)° 



where in the above we used the equahty 11^=1 
us now consider the sum 



. Let 



(21) 



^ a(^) " = 5^ 1 [a(5) = +1] " - J] 1 [a{S) = -1] " 

s s s 

= J2 (1+ [S] + 1- [S] + lo [S]) (1 HS) = +lf - 1K5) = -1] 

s 

= J2 (^lo [S]^f + 1- [S] 1 [a{S) = +l] "-1+ [S] 1 HS) = -1] = 0. 

Adding the above representation of zero to the terms inside the curly brackets in the 
equation fl20l) gives 
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A(m) 




- 1 



U[S]l[aiS) = -l] 











1 — m 


2 


1 - 


1 + m 








v 







1_ [5]l[a(5) = +l] 



(22) 



which is clearly A(m) > for mG [0, 1) and A(m) < for mG (—1, 0]. 

One can show that the function fyim), which we used in the bounding procedure (120 



f l22|) . has the following properties: [3 Tj (i) m > fxi^) when m G (0, 1) and /x(^) > fn when 
mG(-l,0) (/-^(0) = 0) for tanh/3<6(A;); (ii) for tanh/3>6(A;) there exists m* G [-1, 1] \ {0} 
such that fx{m*) = m* . | 



The consequence of Proposition [rvTT] is that the ordered ferromagnetic (FM) phase m^Q 
is a fixed point of f[T7l) (if at all) on/i/ for values of /3 and k which satisfy tanh/3 > h{k). 
This situation leads to the PM/FM phase boundary, depicted in the phase diagram (Figure 
[3]), which for k oo approaches p = 1/2 as 1/2 — p{k) = 0{l/\/k); this is follows from 
the Stirling's approximation of b{k). BN constructed from a single Boolean function x{S) 
saturates this boundary. We note that a similar result, for odd k only, have been conjectured 
using the annealed approximation and multiplexing techniques j26| . 

In the original RBN model, which we have considered in the section llV A 21 the stationary 
state m = 0, g 7^ is for any (3 G (0, oo). Here we explore a possibility for the system ([1]) to 
have the disordered PM (m = 0, q = 0) and two ordered FM {rriy^O, q^O) and SG (m = 0. 



qj^O) states. For lim(_^oo '"^(^) =''^, q = is a fixed point of f|T8|) iff ot{S)}'^ =0 which 
occurs only for balanced Boolean functions, with an equal number of dd in the output. 



Proposition IV. 2 For m = the point q = is a unique stable solution of IfT^) when 
tanh'^ 13 <b{k). 

Proof In order to show this we first define the function T{C) = 



tanh2(/3)E.5n-=i 



sgn[5' • S] which is related to the function via the 
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FIG. 3: Phase diagram for the recurrent Boolean network governed by ([T]). The model is param- 
agnetic (PM) for any distribution of Boolean functions for the noise parameter p > {1 — b{k))/2. 
It can be ferromagnetic (FM) below the boundary (1 — h{k))/2 and there is a possibility of the 



spin-glass (SG) phase when p < {1 — y^b{k))/2. 

equality T(C) = tanh(/3)/^(C). This property follows from the calculation 

k 



T{C) = tanh2(/3)^J] 

k 

= tanh^(/3)^n 

s j=i 



1 + SjSjC 



Sfi 



sgn[5'.5'] 

k 

sgn[^^,]=tanh(/3)/^(C). 



(23) 



Next we define the function g{C) = tanh^(/3) ^^j^n^.^-,^ ^ ^ a{S)a{S), where a is 
an arbitrary balanced Boolean function, and compute the difference A(C) = (T(C) — 
(7(C))/4tanh^(/3) using the same steps as described in the equations (|2n]) - (l22]) . The re- 
sult of this computation is that A(C) > and A(C) < on the intervals C G [0, 1) and 
Cg (-1, 0], respectively, from which the bounds T(C) > F„(0, 0, C) and T{C) < F„(0, 0, C) 
on the same intervals follow. The behavior of T[C] with respect to the inverse temperature 
/3 is the same as of f^, which we described in the proof of Proposition IIV.II but with the 
tanh(/3) being replaced by the tanh^(/3). | 

From the Proposition II V. 2] the case of m = 0, qy^O and finite /3 occurs only (if at all) when 
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tanh^ /3 > b{k). The resulting FM/SG phase boundary (see Figure [3j) approaches p = 1/2 
as 1/2 — p{k) = 0{k~^^^) when A; — )• oo. The a-averages in equations f|T71) - f|T8l) can be 
computed for a uniform distribution over all balanced Boolean functions to obtain m(t) = 
for all t>0, which implies g = tanh^(/3) (^{^)''{^ + 2^)^w^^- "^^^ latter has only one q = 
trivial solution for any finite (3 and develops a second q=l solution only for /3— j-oo. Thus, 
only the model ([T]) with non-uniform distributions over the balanced Boolean functions can 
have the critical behavior as in Figure |3l 

It is interesting that the upper bound b{k) computed here for k odd is identical to the one 



computed for noisy Boolean formulas 2J] . A noisy Boolean formula is a tree in which leaves 
are either Boolean constants or references to arguments, internal nodes are noisy Boolean 
functions (for each function-input there is an error probability p which inverts the function- 
output) and the root corresponds to the formula output. The MAJ-Zc function, which plays 
a prominent role in the area of fault-tolerant computation (FTC) as it allows to correct the 
errors by its majority- vote function {2^, saturates the bound b{k). The idea to have two 
copies of the same system, used in this work only to study initial-states perturbations, is also 



useful for FTC as it allows to compare the noisy system against its noiseless counterpart |25 |. 

The connection of our work with FTC stems from the fact that each site i at time t in 
our model can be associated with the output Si{t) of a k-aij Boolean formula of depth t 
which computes a function of the associated initial states (a subset of {^^(O)}) |ll|- In the 
presence of noise, a formula of considerable depth (large t) loses all input information for 



tanh/3<6(/c) and odd k [2^. This suggests that the upper bound b{k), for odd k, is more 
general and is valid for transitions at all m values identifying the point where stationary 
states depend on the initial states and ergodicity breaks. For k even such general threshold 
is not yet known. 

B. Boolean networks with memory 

In model ([T]) the state of site i at time t depends on its states at previous times only 
indirectly via the sites affected by the state of site i at previous times. These dependencies 
create correlations via the directed loops in the time-space picture of BN (as in Figured]), 
but in the limit of — 00 they become very weak, as was argued in previous works in 
this area ll|]. This allows one to express the observables of interest ([7]) in the closed form 
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(fT7|) -f ll9p . However, in a broad family of models, which includes the Boolean networks with 
reversible computation [28] and gene networks with self-regulation 29|], the state of a site i 
at a time t+1 depends directly on its state at time t. 



1. Random threshold networks 

An exemplar model with strong memory effects, which was used in [30;] to construct a 
model of cell-cycle regulatory network (A^ = ll) of budding yeast, is of the form 

Si{t+l) = sgn[h,{t)-2h] + Si{t)6[h,{ty,2h], (24) 

where ^(t) = Ej=i (1 + -^i^ (t)) and G {-1, 1}. Mean-field theory (iV — 7- oo) was de- 
rived 1J| using the annealed approximation in a variant of this model, where the interac- 



tions {^j} were randomly distributed P{^j = ±1) = 1/2. Significant discrepancies between 



the theory and simulation results has been pointed out for integer h values (in this 
case it is possible that 2h = hi{t)), which was attributed to the presence of strong memory 
effects. Refinements of the annealed approximation method improved the results obtained 



only slightly SjJ, |32| but break down in most of the parameter space. 



This model f l24j) can be easily incorporated into our theoretical framework. The result of 
the GFA fllSp for this process (with thermal noise) can be obtained by replacing the average 



(■ ■ ■ ) by (■ ■ ■ ) and the probability function Pa{S{t+l)\Si{t), .., Sk{t)) by 

/3S(t+l ) {sgn [/i (t)-2/i}f S(t) 5 [/i (t ) ; 2/i] } 

where A(() = E'.i&(l+S,(()). 

In the case of /igM, the probability function fl25|) is independent of S{t) and equations 



T7|) - (|T9|) have the same structure as model (12^ : the a-averages a{S) and a{S)a{S) 



are replaced by the averages sgn[h{t)—2h] and sgn[hit) — 2h] sgB.[h(t)—2h] , respectively. 
The equation for m{t) recovers the annealed approximation result 14] (using the relation 
b{t) = {l+m{t)) / 2) . In Fig. |l|^a,b), we plot our analytical predictions for the evolution of m{t) 
and C(t + t^,t^) against the results of Monte Carlo (MC) simulation which use flM|) . The 
correlation function C(t + tiu,tiu), in the limit of t — )■ oo, — )■ oo, approaches the stationary 
solution of the overlap function f|T9l) with increasing as predicted (Fig. Hll^b)). 

The situation is very different when h E 1^. The magnetization m{t) = P{S)S{t), 
where P{S) is a marginal of (fT^ with P^-^P^, is no longer closed as in (ITTj) . but depends 
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FIG. 4: Evolution of the magnetization (m = m[t)) and correlation (C = C{t+tyj,tyj)) functions 
with time t is governed by ()24p . Theoretical results (lines) are plotted against the results of MC 
simulations (symbols) with A^ = 10^. Each MC data-point is averaged over 10 runs. Error bars are 
smaller than symbol size. Evolution of m (a) and C (b) for /iGM. In (b) we plot C for /i = 0.5 and 
k = 2>. 




02468 10 2 4 6 8 10 

t t 



FIG. 5: Evolution of m = m{t) (a) and C = C{t+tyj,tyj) (b) for /iGZ. In (b) we plot C for /i = 
and k = 2. Theoretical results (lines) are plotted against the results of MC simulations (symbols) 
with N=W^. Each MC data-point is averaged over 10 runs. Error bars are smaller than symbol 
size. 

on 2*~^— 1 macroscopic observables (all magnetization, all multi-time correlations). Thus the 
number of macroscopic observables that determine the value of m{t), or any other function 
computed from f|T5|) . grows exponentially with time. Annealed approximation results [ij] for 
this model when /i G Z are only exact up to t < 2 time steps (the equation for 6(1) = (l-Hm(l))/2 
in our approach and in [14J are identical) and deviate significantly from the exact solution 
at later times (Fig. O^c)). A typical evolution of the correlation function C{t+ty^,tw) in the 
system ( I24p when h ^ Z is shown in Fig. El^d). 
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V. DISCUSSION 



We applied the generating functional method to analyze the dynamics of recurrent 
Boolean networks. The analysis resulted in a coupled set of recursive equations for a small 
number of macroscopic observables that provide an exact description of the dynamics in a 
broad range of Boolean networks. Based on the analysis we also showed that for a large class 
of models the annealed approximation does provide exact results for both magnetization and 
overlap order parameters; although results for correlation between states at different times 
are generally incorrect. However, it turns out that in models where the state of spin at 
time t depends on its state at previous times directly the annealed approximation always 
fails. This is due to the fact that approximation does not take into account the correlations 
which are very strong in such models. Comparing the two transition probabilities ([2]) and 
fl25p for the processes without and with memory, respectively, it is clear from the single-spin 
trajectory equation f[T^ . that as soon as there is an explicit dependence of a spin on its 
states at previous times, an exponential (in time) number of macroscopic observables will 
be required. Furthermore, also in models where this approximation works well it is useful to 
know the two-time correlations as it provides us with insight into properties of the stationary 
states. 

When one considers systems with thermal noise, the suggested framework provides addi- 
tional new and interesting results. We have computed the noise threshold above which the 
system is always ergodic and critical noise levels where phase transitions occur. Here the 
two-time correlation function is particularly useful as it allows one to compute the Edwards- 
Anderson order parameter q, used to detect the spin-glass phase. 

One of the remaining questions is to provide an example of a model (or show that it does 
not exists) with a phase diagram as in Figure O The direct computation of q for all balanced 
Boolean functions with k > 3 inputs is possible for small k but soon becomes intractable 
as the number of such functions grows exponentially with k. The other important question 
is to find a systematic way to generate good approximations for the dynamics with strong 
memory effects. Although the theory developed in this work can describe the dynamics of 
such models exactly it can be used only for relatively short times due to the rapid increase 
in the number of macroscopic order parameters. However, the equations of our theory can 
be used to check the quality of approximations used in such models in all future works 
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and possibly serve as a starting point for such studies. The work undertaken here can be 
extended in a numerous ways. For instance, one can easily adapt the framework developed 
here to study Boolean networks with inhomogeneous connectivities 



different noise models 



33| and to examine 
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Appendix A: Computation of the disorder averages 



In this section, we outline the calculation which takes us from the definition of generating 
functional (E]) to the saddle-point integral ([8]). The starting point of this calculation is the 
generating functional 



(Al) 



{Sl{t),S^{t)} 



X exp 



fc=0 ^1 0=1 



where in the above we have defined the field- variables h]{S^{t)) = 
J2h jfe ^ii,...,ifc'^«('5'i'| (''');•••; Removing these fields from equation ( 1A1I) via 
the integral representations of unity 



n nn 

fc=0 i=l 7=1 



d/i7(t) d/j.7(t) i/,7(t)[/^7(t)_fe7(5^(t))] 

277 



(A2) 



gives 



itnax ^ 2 

5^ P(5i(0), 5^(0)) exp -i^EE^'W^'W 

{S}(t),S^{t)} L fcO ^1 T=l 

2cosh[/3/iJ(t)] 



(A3) 



27r 



X nnn 

t=0 i=l 7=1 

i"m,ax 1 2 

t=0 i=l 7=1 



a,(57^{t),...,S7Jt)) 



(A4) 
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Now we can average out the disorder in (lA4p 



tmax-l N 2 N 



n nn n ^"^^^^^ .^^.{s:;^,...,^;;^) 

t=0 i=l 7=1 ii 



^i=l 1C[N] Al 
N 

X j]PK) n 

oii il,---,ik 

^i=lIC[Af] Al 
N _ 

X n e 

ii,...,ik 



}' 


1; Y: Ar 




l'C[N] 







1; Y 






I'C[Af] 



X exp 



1 



dWi : 
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and use this result in the generating functional flASP to obtain 



(A5) 



{Sl{t),Sfit)} 



tjnax ^ 2 

fc=o i=i r=i 



t=0 i=l 7=1 

AT 



X 



2cosh[/3h]{t)] 

^^^/^e-}-p[^'/w)}/_>^ 

-5^ 5(0;-.;.) n W)-h.(t)) 
i=l I t=o J 



N 

i=l V, 

-I [ i"max 1 k 



(A6) 



where in the above we have defined the vectors hj(t) = {h}{t), hf{t)) and Si{t) 
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{Sl{t), Sf{t)). Using the unity representations 



(A7) 



gives 



X 



expiv[i J2 Pi{Sm)Pi{Sit)}) + i /"{dh(t)}dc^^]({h(t)},a;)^]({h(t)},u;) 



{k 
llP{{SAt)}) 
j=i 



(A8) 



i?=0 ^1 1=1 



X J] P(5i(0),52(0))exp 

27r 2cosh[/3/iJ(t)] 



X nnn 

t=0 1=1 7=1 
N 



i=l 



2tx 



Equation flASp gives one the saddle-point integral ([S]) if write its site-dependent part in the 
exponential form 

N 



exp 



Elog /{dhKt)dh,Kt)} r^M[{5,(t),h,(t)}|{h,(t)},a;,,R: 

i=l i^.itw 



m 



where the definition 



M[{5.(t),h,(t)}|{h,(t)},a;„{^7(t)}] 



P{Sl{0),Sfmew 



t=0 7=1 

xe 



ihj(t)hj{t) e ' 



,/357(t+l)fe7{t) 



2cosh[/3/i7(t)] 

Q({h,(i)},cj,)+i'^>-if'({'S',(i)}) 



(A9) 



is used with the shorthand / {dhi(t) dhi(t)} = HtlT 117=1 / 



ima:.-lT-r2 T d/r7 (t) dfe7(t) 



27r 
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Appendix B: Solution of the saddle-point problem 

In this section, we show that the conjugate order-parameter function P{{S{t)}), which 
is governed by the equation (fT2i) . is a constant function. In order to do this, we first rewrite 
the (disorder-averaged) path-probability ([3]) as follows 



p[{s\t)y,{s\t)}] = ^p{s\o),s\o)) 



^ { imax 1 2 



""Eli G 32cosh[/3EL...^l,....."^(^nW'---'^i.W)] 



x5 



I'c[Af] 

N 



1=1 



{dh,(t) dh,(t)} 



27r 



X • 



^max 1 2 

nn 

<=0 7=1 



0/357 (t+i)h7{t) 



/ {dh.(t) dh.(t)} r ^ Ha{5,(t),h,(t)}|{h,(t)} 



2cosh[/3/iJ(t)] 



(Bl) 



p{sm,sm) (B2) 



where in the above we have averaged out the connectivity disorder only, i.e. ( 
^{A\} riicfAf] + (l~]^)'^Aj;o| (' ' ' ) ^^Le definition of single-site measure Sj i 

clear from equation (IB2p . 
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In the next step, one notes that the effect of the Fourier transform 



{dh,(t) dh,(t)}/ S,[{5,(t),h,(t)}|{h,(t)},a;, 



xe 



J {dh,(t) dh,(t)} ^ e-'(°-^n,.....^n,.....)p(S,l(0),5f(0)) 



X 



X 



tmax 1 2 

nn^ 

<=0 7=1 
p/357 (4+1)^7 (t) 



2cosh[/3/i7(t)] 



a.(57^{t),...,S7jt))] 



'imai;-! 2 ^13 1+1)91 (t) 

11 ll2cosh[/3^7(t)] 



AT 



0? 



P(5,i(0), 5,2(0)), 



(B3) 



on the function Sj is to replace the Boolean function on site i with a constant function 
9^{t) G { — 1,1}. With this in mind we can define the average site-perturbed path-probability 



(B4) 



i=l 



N 



X 



N f N 

E n /{dh,(t) dh,(t)} 

i=i y -J 

^H,[{5,(t),h,(t)}|{h,(t)},a;, 



X /{dh,(t) dh,(t)}/ ^S,;[{5,(t),h,(t)}|{h,(t)},u;, 



xe 



To show that the object in flB4p indeed defines a probability measure we note that it is 
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clearly a positive semi-definite and the normalization of (lB4p can be established as follows 

(B5) 



1 ^ 



{Si(i),S2(t)} i=l 



E n E /{dh,(t)dh,(t)} 



^ dw, 



^S,[{5,(t),h,(t)}|{h,(t)},u;, 



X 



J2 JidHt) dHmj^^ S,[{5.(t),h,(t)}|{h,(t)},a;, 



{Si(t)} 

AT AT 

N 



^ N N 



1=1 

AT AT 





-{A^J— 




i;E 




o;E^i 


I'C[Ar] 




IC[N] 



»=1 j^i {Aj}IC[Af] 



l;E^r 

I'C[Af] 



IC[Ar] 



exp 



o;E^i 

IC[Af] 



1 



(1- 



\Ar*-^l 



Af-l 



(1- 



-AT 



Since for — t- oo the inverse of Za in flB5D grows as e , we conclude that 



limAT^^oo Z]{si(t),s2(t)} ^ Zlili {S'^(^)}]u^^e = 1- Alternatively, the calculation 
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in (IBSp can be carried out differently using results from the appendix \M 



1 ^ 



{Si(t),S2(t)} i=l 
y-1 N ( N 



E n E /{dh,(t) dh,(t)} 

r duj- {^^.a 

X / ^S,[{5,(t),h,(t)}|{h,(t)},u;J ^ 

^ |{dh.(t) dMt)}J^^ E,[{s,it)Mt)}\{Ht)},. 



{S.(i)} 

xe" 



iE«""~'h,(t).e(t)-iL^, 



N ^ /{dP'dP'dl^'dl]'} e^*[{-P''^''^'.f^'}] \ / 



where in the above Mj-average is generated by the site-dependent ver- 
sion of ([9]). Computing the integrals in (\B6\i by the saddle-point 



method we find that limTv^oo E{siw,s2(t)} ^ Eii ^[{-^^ W); {-^^Wllu^^e 
/{dh(t)} dui ^l{{h{t)} , Lo) e'^^'S)""' h{t).e(t)-iuj ^ where the order-parameter function Q. 
is defined in flTS]) . but according to the calculation in f lBSp this also equals unity. Thus, by 
using this result in the saddle-point equation ( TT^ . we find that the equality P{{S{t)}) = ik 
holds. 
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